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ABSTRACT 

Galaxy spectra are a rich source of kinematical information since the shapes of the 
absorption lines reflect the movement of stars along the line-of-sight. We present a 
technique to directly build a dynamical model for a galaxy by fitting model spectra, 
calculated from a dynamical model, to the observed galaxy spectra. Using synthetic 
spectra from a known galaxy model we demonstrate that this technique indeed recovers 
the essential dynamical characteristics of the galaxy model. Moreover, the method 
allows a statistically meaningful error analysis on the resulting dynamical quantities. 

Key words: Galaxies: kinematics and dynamics-Gelestial mechanics, stellar dynam¬ 
ics 


1 INTRODUCTION 


The construction of a dynamical model for a galaxy de¬ 
pends on two kinds of data. The photometric measure¬ 
ments yield, after deprojection, the mass density. Assuming 
a mass-luminosity ratio, the Poisson equation then leads to 
the gravitational potential. The kinematics are obtained by 
analysing spectral information, with the aim to characterise 
the line-of-sight velocity distribution (hereafter LOSVD) 
(j){vp,Xp,yp), which is the probability of finding a star with 
projected velocity Vp at a position {xp,yp) on the sky. 

When the signal to noise ratio {S/N) of the data is 
rather low, a LOSVD can be parametrised as a simple 
Gaussian. There exist a host of methods to determine the 
best fitting Gaussian LOSVD, such as the Fourier quotient 
meth od (S imkin 1974, Sargent et al 1977), modified by Win- 
sail (1991) to incorporate also non-Gaussian LOSVDs, the 


cross correlation method (Tonry & Davis 1979), the Fourier 
fitting method (Franx, Illingworth & Heckman 1989) and 
the Fourier correlation quotient method (Bender 1990). 

Motivated by the increasing quality of the data, new 
methods have been developed recently that are able to quan¬ 
tify symmetric and antisymmetric deviations of the LOSVD 
shape from a Gaussian (e.g. van der Marel & Franx 1993 
and Gerhard 1993). A LOSVD is then parametrised by the 
average projected velocity (vp), the projected dispersion Up 
and the lowest order coefficients of a Gauss-Hermite expan¬ 
sion of the LOSVD, hs and h^. These can subsequently be 
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used in a fitting routine to construct a dynamical model that 
best represents the kinematical data. 

Another path leading to estimates of the shape of 
LOSVDs employs non-parametric fitting procedures such as 
the Wiener filtering method (Rix & White 1992), the max¬ 
imum entropy method (Statler 1991), the unresolved Gaus¬ 
sian method (Kuijken & Merrifield 1993) and the Bayesian 
fitting method (Saha & Williams 1994). These give a totally 
unbiased estimate of the LOSVDs. 

It is clear that all these modeling techniques require 
two fits. One fit is needed to extract the kinematics from 
the spectra, and a second fit then produces a dynamical 
model. In this paper we propose a method to construct a 
distribution function that operates directly on the available 
spectra, after they have undergone all the necessary steps of 
a standard data reduction. This probl em is well defined, as 
is shown by Dejonghe & Merritt ( 1992 ). The essentials of the 
method are explained in the next section. In section 3, we 
test the method on a set of spherical anisotropic Plummer 
models. 


2 THE MODELING PROCEDURE 

A galaxy spectrum, obtained at a given position {xp,yp) on 
the sky, is the integrated light of all the stars along the line- 
of-sight. The spectrum of each contributing star is Doppler 
shifted according to its velocity relative to the observer. If we 
assume that the bulk of the light of the galaxy is produced 
by stars of approximately the same spectral class (in the 
case of an elliptical galaxy, most of the light originates from 
the G to M giants) we can write the following form for a 
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galaxy spectrum at {xp,yp) 


g{xp,yp,X) 


j 0(wp, Xp,yp) s[X{vp)) dvp 


( 1 ) 


which expresses the spectrum as a sum of products of prob¬ 
abilities : the probability of finding a star with projected 
velocity Vp at a position {xp,yp) and the probability that 
such a (standard) star emits a photon with Doppler shifted 
wavelength A, according to its spectrum in the rest frame 
s(A). We denote the rest wavelength by X{vp). For velocities 
which are small compared to the velocity of light one has 
the Doppler shift formula 


Xp X A(up) __^ A 

c X{vp) X{vp)' 


( 2 ) 


Obviously, (|^) can be written in the form of a convolution us¬ 
ing logarithmically rebinned spectra, which also eliminates 
the choice of a dimension for the wavelengths. Hence 


g{xp,yp,\nX) = J (j>{vp,Xp,yp) s{lnX - ^)dvp. (3) 

It is always possible to write the distribution function 
(hereafter DF) F’(I), where I are the isolating integrals of 
the motion, as an infinite sum of basis functions ^*(1) with 
coefficients d 




( 4 ) 


In the sequel, we will consistently use a summation conven¬ 
tion. These basis functions are the DFs of simple dynamical 
models, called the components, that form a complete set. 
Denoting by Vx and Vy the velocity components in the plane 
of the sky, and by 2 the coordinate along the line-of-sight, 
we find formally 


/ OO 

•OO 

Obviously, the LOSVD is linear in the DF, and can also 
be written as a weighed sum of component LOSVDs 
(jf{vp,Xp,yp) with the same coefficients a as the DF. Like¬ 
wise, we obtain component spectra 

g\xp,yp,lnX) = J (j>''{vp,Xp,yp) s(lnX - ^) dvp. (6) 

We label with the composite index n the fe-th pixel of the 
l-th galaxy spectrum (we label each spectrum and its corre¬ 
sponding line-of-sight with a number 1) and denote its value 
by gr„. A similar notation applies for the spectra of the i-th 
component, g!^. It is easy to see that ideally 

gn = Cigl, n = l,...,N, (7) 


dz F{T)dVxdVy. 


( 5 ) 


with N the total number of pixels considered. 

Since we can only handle a truncated expansion of the 
DF, we will have to use a limited number m of compo¬ 
nents and try to find those coefficients d that give rise to a 
model that best represents the data. This we do by using a 
Quadratic Programming minimisation technique (Dejonghe 
1989). The stellar template spectrum is normalised to one 

/ OO 

s(ln A) din A = 1, (8) 

•OO 

where the function s(ln A) is set to zero in the region outside 
the data. The aim is to model the observed spectra with (j^. 


Since the component LOSVDs 4’''{xp,Xp,yp) are normalised 
to their projected mass density Pp{xp,yp), the calculated 
component spectra g''{xp, yp, In A) will also be normalised to 
Pp{xp,yp). In order to make (^ correct, in the sense that 

P{xp,yp) = dpl{xp,yp), (9) 

the observed galaxy spectra g{xp,yp,\nX) must be nor¬ 
malised to the projected mass density Pp{xp,yp). However, 
an observed galaxy spectrum cannot be written exactly as 
the convolution of a template spectrum and a LOSVD since 
such a convolution has two rapidly decaying wings in both 
regions where the template spectrum is unknown and by de¬ 
fault set to zero. Observed galaxy spectra are also unknown 
and set to zero in the same regions as the template spec¬ 
tra. So by normalising the galaxy spectra to the projected 
mass density a small error is introduced, because the two 
small but unknown wings, required to make the convolution 
formally exact, are neglected. 

Usually the continuum of a galaxy spectrum and the 
continuum of a weighed sum of component spectra will not 
be the same. That is why the continua of the galaxy spectra 
and the stellar template spectrum have to be subtracted be¬ 
fore they are normalised. For simplicity, we still use the no¬ 
tation g{xp, yp, In A), now for a continuum subtracted galaxy 
spectrum, and s(lnA) for the continuum subtracted stellar 
spectrum. 

The quantity that has to be minimised is 


= w" (gri - d g^Y , (10) 

with to" = 1 /ct^ and an the noise on gr„, which is well mod¬ 
eled by Poisson noise. This expression can be expanded into 
the form 

= Ci D’'Rj - 2p'd + e (11) 

or, in matrix notation 

= c*Dc — 2p*c-I-e. (12) 

The matrix D with elements 

= w-gi,gi (13) 


is called the Hessian matrix and the vector p has components 


P = W Pngn 


(14) 


The scalar constant e is given by 

e = w'^gngn- (15) 

This x^ is quadratic in the coefficients. It has to be min¬ 
imised under the linear constraint that the DF has to be 
positive 

dF\l) > 0. (16) 


The positivity is tested on a grid in phase space. This is a 
typical problem of Quadratic Programming. Moreover, the 
X^ has the statistical meaning of a goodness of fit so its 
expectation value equals 

ix^) = N -m-1. (17) 

Other observational data can be added to this . However, 
it looses its meaning as a goodness of fit when other weights 
than the errors are used in In appendix C we present 
the error analysis on the obtained coefficients. 
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Figure 1. The LOSVDs of a radial Plummer model with q = 1 (top) and a tangential Plummer model with q = —6 (bottom), both 
with M = 5 X 10^^M© and rc = c = 5kpc, for rp = 0, 1, 5 and lOkpc (from left to right). 


3 PUTTING THE METHOD TO THE TEST 

To test our modeling procedure we make use of synthetic 
galaxy spectra with simulated Poisson noise. These are cre¬ 
ated with spherical anisotropic Plummer models (Dejonghe 
1987). The synthetic spectra are then used as input for 
our modeling procedure. As components we use generalised 
Fricke models. 


3.1 The Plummer model 


The Plummer potential can be written as 


r/)(r) 


GM 1 

a/ 1 + (r/rc)2 


(18) 


with Vc the s o-cal led core radius and M the total mass. We 
take a/ tpiO) = a/ GM/tc as the unit of velocity so the cen¬ 
tral value of the potential becomes unity. The density of the 
anisotropic Plummer models can be written as 

p{r) + (19) 


with q and c real numbers. This function is an augmented 
mass density that is instrumental in the construction of the 
2-integral DF F{E,L) (Dejonghe 1986). Inserting the func¬ 
tional form for the Plummer potential in this expression, 
with Vc = c, one sees that all Plummer models have the 
same density distribution 

p(r) = (1 + (r/cf)-^!\ ( 20 ) 


We will approximate the anisotropic Plummer model with 
Fricke components with augmented mass density 

p{r) ='i}(rY {r! cf^ ( 21 ) 


with a and /? real numbers. The characteristics of both the 
Plummer models and the Fricke components can be calcu¬ 
lated using a more general expression for the mass density 


/' 

2/3 / /, 

,A 2\ 


p{r) = y{ry / 

0 '+(: 

/) 

(22) 


with a, P and 7 real numbers that satisfy 


a > 3 - 27 (23) 

to keep the total mass finite. Choosing /3 = 0, a = 5 — 
q and 'y = q/2 one recovers the density of the Plummer 
models. The choice 7 = —/3 yields the density of the Fricke 
components. The advantage of the above-mentioned family 
lies in the fact that most of the kinematics can be calculated 
analytically. This is for instance the case for the anisotropic 
velocity moments (see Appendix A) 


P2n,2mir) = C{m,n)^{ry 


(r/c) 


2/3 


(1 -b (r/c)2)/3+3' 

(r/cf , 


with 

C{m,n) = 


2Fi{-m,-y + p-p+l, 

2™+" r(n -b l/2)r(Q; -b l)r(/3 + m+l) 
ypr r(/l -b l)F(fr -b m -b n -b 1) 


(24) 


(25) 


This general expression yields the following expressions for 
the radial and tangential velocity dispersions, 


2/3 


p{r)aUr) = P2,o(r) = Y^^(r)^+“ 

p{r)cjl{r) = p{r)al{r) = \pq, 2 {r) = x 

2 1 + a 


(r/c) 


2/3 


(1 -b (r/c)2)/3+7 


1 - 


7 -b /? {r/cf 
1 -b/3 (1 -b (r/c) 2 ) 


Binney’s anisotropy parameter then reads 
_ -P + '^{r/cf 


P{r) = 1 - -I = 


1 -b {r/cY 


(26) 


(27) 


(28) 


For a detailed derivation of the anisotropic velocity mo¬ 
ments, we refer to appendix A. 

One easily sees that the Fricke components have a con¬ 
stant anisotropy 


P{r) = -p. 


(29) 


The tangentially anisotropic components have zero central 
density whereas the radially anisotropic ones have a central 
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density cusp. The isotropic components have a nonzero and 
finite central density. For the Plummer models 


P{r) 


q jr/cf 

2 1 + (r/c)2 ■ 


(30) 


We see that /3(r) and q have the same sign and since —cx 3 < 
/3(r) < 1 the allowed range for q is —(Xj < q < 2. The true 
nature of the orbital structure of a model becomes clearly 
visible in the kinematics at radii r > c. 

A LOSVD of a tangential Plummer model along a line- 
of-sight close to the center is narrow and peaked. Mov¬ 
ing the line-of-sight outward, the LOSVDs become progres¬ 
sively broader and ultimately obtain a bimodal shape. The 
LOSVDs of a radial model, on the contrary, will be broader 
towards the center and narrower at large projected radii. 
As an illustration. Fig. ^ shows the LOSVDs of a radial 
Plummer model with q = 1 (top) and a tangential Plummer 
model with q = —6 (bottom) for Vp = 0, 1, 5 and lOkpc. 


3.2 The synthetic spectra 

We want to construct spectra with (|^ so it is important 
that we can calculate the necessary LOSVDs efficiently. It 
turns out that for a subset of the models one can perform 
the integrations of the DF over the velocity components in 
the plane of the sky analytically. This yields lp{r, rp,Vp), the 
probability of finding a star with projected velocity Vp at a 
distance r from the center of the cluster and on a line-of- 
sight at a distance Vp in the plane of the sky from the center. 
The following expression is valid only when /3 is a positive 
integer and /3 + j a, negative integer (see Appendix B) : 


Mrr V) = 1 r(a + l) (r/c)^^ 

r(a + 1/2) (1 + (r/c)2)/3+7 ^ 

-(/3+7) • 


Ed). 


(-/3)i 

+3 (*+/)! 


2 \ (-(i+j) + a-i-/2) 

I _ 


2Tp{r) 


Here 


(sin 77 ) 


1 vt 




{r/cf 


1 -I- (r/c)2 


(31) 


(32) 


and T](r, Tp) is the angle between the spherical radial di¬ 
rection and the line-of-sight at a certain point of the line- 
of-sight. The integration along the line-of-sight has to be 
performed numerically. The LOSVDs of the models we use 
are particularly simple. In the case of Fricke components the 
above expression is reduced to a sum over i, for Plummer 
models a sum over j remains. LOSVDs of models that do 
not have positive integer or negative integer -I- 7 have to 
be calculated numerically using (^. 

We calculate spectra out to about two effective radii. 
For a de Vaucouleurs law, a distance of two effective radii 
encloses approximately 69% of the total luminosity. For a 
Plummer galaxy this corresponds to rp = 1.49c. 

Besides a LOSVD, one also needs a stellar template 
spectrum to construct a galaxy spectrum. We use the spec¬ 


trum of the K giant HR2425 in a wavelength range from 
3900A to 5400A. This spectrum is rebinned logarithmically 
so that one pixel corresponds to a velocity interval of 52 
km/s. The Gaussian width (uinstr) of a line that is intrin¬ 
sically a delta function is about 150 km/s. The synthetic 
spectra are calculated with (^ and noise is added according 
to a specified S/N ratio. 


3.3 Testing the method 

To test if our modeling procedure can successfully perform 
the inversion of the integral equation (^ we offer the pro¬ 
gram noisy synthetic spectra of a Plummer galaxy. We as¬ 
sume that the Plummer galaxy contains only stars of the 
same spectral class. In that case it is not necessary to sub¬ 
tract the continuum of the spectra. We also use the stellar 
spectrum with which the synthetic spectra are created as a 
template spectrum. These assumptions are needed to avoid 
problems with bad continuum subtraction, template mis¬ 
match etc. that have only to do with the data reduction and 
are common to all methods. Due to the normalisation of the 
template spectrum and the LOSVDs, the synthetic spectra 
will have automatically the proper normalisation (see the 
discussion about the normalisation of spectra in the para¬ 
graph below equation (^)). 


3.3.1 A tangential q = —2 Plummer galaxy 

We first consider a mildly tangential q = —2 Plummer 
model. The model has a core radius = c = 6.75kpc and 
a total mass of 6.5 x10^^Mq. This galaxy is placed at a dis¬ 
tance of 206265 kpc so that one arc-second (”) on the sky 
corresponds to a distance of 1 kpc. Two effective radii cor¬ 
respond to approximately 10” so we generate spectra for 
Vp = 0 ”, 0 . 5 ”, 1 ”, ... , 10 ”. We only use a wavelength inter¬ 
val of about 170A containing the Mgl lines around 5180A 
and a strong Fell fine at 5270A. During the fit, the correct 
mass and Plummer potential are adopted. 

The augmented mass density of this Plummer model 
can be written as a binomial expansion, which shows that it 
can in principle be fitted exactly in terms of Fricke com¬ 
ponents with a = 7. We first offer the routine noiseless 
data. A library is set up that contains the Fricke compo¬ 
nents with a = 4, 5, 6, 7, 8, 10 and 12 and with positive 
integer /3. These components are chosen because they span 
a wide range in dynamical behavior. The models with high 
Q-values are more centrally concentrated. The higher the 
value of P the more tangentially anisotropic the component 
becomes. When dealing with real data, one has to make a 
reasonable choice of components. In the following we will 
only use Fricke components with positive integer /3. This is 
just to save calculation-time ; the LOSVDs of these compo¬ 
nents can be calculated analytically. 

The ’redundant’ components do not contribute signif¬ 
icantly to the result (they are not chosen or else given a 
very small coefficient) and the correct components have the 
right coefficients. Thus, the input model can be successfully 
recovered, which is of course a prerequisite. 

We now consider data with a signal to noise ratio of 
S/N ~ 80 in the center, dropping to S/N ~ 25 at two ef¬ 
fective radii. Here and in the following, all signal to noise 
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ln(X(A)) v(km/s) 


Figure 2. The synthetic spectra of a tangential q = —2 Plummer galaxy (black) with central S/N ^ 80 for four different lines-of-sight. 
The scaling of the ordinate is arbitrary. Over-plotted in grey is the best fitting spectrum (the Fricke components with a = 4, 5, 6, 8, 9, 10 
and 12 are used). To the right of each spectrum, we plotted (0(theo)-0(fit))/^(theo)7n.aa; with 0(theo) the Plummer LOSVD, ^(fit) the 
fitted LOSVD and ^(theo)-max the top value of the Plummer LOSVD. One can see that the deviations between the Plummer LOSVD 
and the fitted LOSVD are nowhere larger than 1 to 2 % of the top value. 





Figure 3. The kinematics of a tangential q = —2 Plummer galaxy (black). Clockwise are plotted the projected velocity dispersion 
C7p(rp), the projected mass density Pp{rp) (in dimensionless units), the tangential velocity dispersion o-(j){r) = o-^^r) and the radial 
velocity dispersion (Jr(r). Shown in grey are the kinematics derived from the best fitting model (using data with central S/N ^ 80 and 
Fricke components with o; = 4, 5, 6, 8, 9, 10 and 12). The deviations are of the order of 5 km/s. The errorbars give a fair idea of the 
actual errors on these derived quantities. There is one active constraint. 
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Figure 4. The DF of a g = —2 Plummer galaxy (black) and the fitted DF (grey), with errorbars, shown for L = 0 (left) and L = Lmax{E) 
(right). The fitted DF results from a fit to data with maximum S/N fa 80. Fricke components with « = 4, 5, 6, 8, 9, 10 and 12 are used. 
There is only one active constraint. 


values are expressed per pixel. The Fricke components with 
a = 4, 5, 6, 7, 8, 10 and 12 again serve as components. It 
appears that the spectra constrain the LOSVDs to a high 
degree. The routine still uses the correct components and as¬ 
signs to them approximately the correct coefficients but now 
also uses a few other components to reach a minimum x^- 
A more objective test is to fit the data without incorporat¬ 
ing the components that can represent the Plummer model 
exactly in the library. The spectra are now fitted with a li¬ 
brary containing the Fricke components with a = 4, 5, 6, 8, 
9, 10 and 12. This also gives an accurate approximation of 
the LOSVDs (see Fig. ^). Although they are not included 
in the fit, the observables o-p{rp) and Pp{rp) are very well 
recovered (see Fig. 

One can also use the errors on the coefficients to esti¬ 
mate the uncertainty on the DF itself. Fig. ^ shows the DF 
of the q = —2 Plummer model and the DF that results of 
the previous fit. The left plot shows a cut through the DF 
at L = 0, the radial orbits. The right one displays a cut 
along L = Lmax{E), the maximum angular momentum an 
orbit with a certain energy can have, which is the locus of 
the circular orbits. It can be seen that the fitted DF closely 
follows the Plummer DF. Only in the center, where the DF 
is rather peaked, is there a noticeable difference. 

When spectra are used with S/N ~ 55 in the central 
regions, dropping to S/N ~ 15 at two effective radii, the 
method still yields a good approximation to the Plummer 
model (see Fig ^ and Fig. 0 compares the Plummer DF 
with the fitted DF. There are no active constraints. The 
deviations between the kinematics of the fitted model and of 
the Plummer galaxy are larger but the errorbars, indicating 
the uncertainty on the coefficients, are accordingly larger. 

The analysis of about 20 spectra, each spanning a spec¬ 
tral region of ~ 200A, with the limited number of compo¬ 
nents used here, takes only a few hours on a Pentium ll (this 


is the time needed to calculate all the component spectra 
and to do the fit). 

3.3.2 Testing the errorbars 

To test if the errorbars give a reliable idea of the 1 u uncer¬ 
tainty on the coefficients, we construct 40 sets of 21 synthetic 
spectra with rp = 0”, 0.5”, 1”, ..., 10” and a S/N « 80 at 
rp=0''. The spectra are generated with the q = —2 Plummer 
model of the previous paragraph. 

To each of these sets a model is fitted, using the Fricke 
components a = 4, 6, 7, 8 and 10 with positive integer (3 
and assuming the correct Plummer potential and mass. We 
then calculate for rp = 0”, 1”, 5” and 10” and for selected 
Vp the mean and the dispersion (this we call the experi¬ 
mental dispersion) of the LOSVDs of those 40 models. The 
average LOSVDs tend towards the true LOSVDs and the 
’experimental’ and the theoretical Icr-errorbars essentially 
coincide, as can be seen in Fig. Hence, the errorbars are 
a reliable measure for the possible spread between the input 
model and a typical fit, as caused by the adopted Poisson 
noise. One also concludes that if the used component library 
is sufficiently complete the systematic errors are negligible, 
compared to the random errors caused by the Poisson noise. 

3.3.3 A tangential q = —6 Plummer galaxy 

This is a fairly tangential model. With the choice rc = c = 5 
kpc, the bimodal structure will reveal itself at rp > 5”. The 
mass is 5 X 10^^M©. Two effective radii now correspond to 
7.5”. Spectra are generated for rp = 0”, 0.5”, 1”, ..., 7.5”, 
once with a central S/N ~ 80 and once with a maximum 
S/N K, 55. This Plummer model can be written exactly as 
a weighed sum of the four Fricke components with a = 11. 

The Fricke models with a = 7, 9, 10, 11, 12, 13 and 
14 serve as components. The routine does not recognise the 
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Figure 5. The same as Fig. g but for a fit to data with central S/N ps 55, using the same Fricke components. There are no active 
constraints. 






Figure 6. The same as Fig. 1^ but for a fit to data with maximum S/N 55, using the same Fricke components. There are no active 
constraints. 
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Figure 7. The same as Fig. Wbut for a fit to data with maximum S/N ^ 55. 


Fricke components with a = 11 as the correct ones. They 
appear in the weighed sum with a coefficient that is totally 
different from the one that would be expected from the bino¬ 
mial expansion of the augmented mass density. It is possible 
to obtain a very good fit to the data (see Fig. ^ and Fig. 
for data with a central S/N ~ 80). The LOSVDs (includ¬ 
ing the bimodal shape at large radii) and the kinematics are 
very well reproduced. Fig. O compares the Plummer DF 
and the fitted DF. There are no active constraints. 


3 . 3.4 ^ radial g = 1 Plummer galaxy 

A radial model has much broader LOSVDs in the center 
than a tangential model with the same mass. In the outer 
regious the inverse is true. As a consequence, the outer spec¬ 
tra of a radial model still contain some details that are to¬ 
tally smoothed out in the tangential models. Therefore the 
inversion of the convolution integral (^ is better behaved 
and the LOSVDs of a radial model are better constrained 
by the data at large rp then those of a tangential one. This 
Plummer model has = c = 5 kpc and a total mass of 
5 X Spectra are generated for rp = 0", 0.5", 1", ..., 

7 . 5 ", once with a central S/N ~ 80 and once with S/N ~ 55. 
As components we use Fricke components with a = 4, 5, 6, 
7 and 8. 

The LOSVDs of this model cannot be written exactly 
as a weighed sum of Fricke components. Nonetheless, the 
LOSVDs are recovered very well in both cases (see Fig. 
and 1^. This is remarkable, because all componeuts are ei¬ 
ther isotropic or tangential. The only way to construct a 
radial model with these components is to give tangential 
components a negative coefficient to suppress the number of 
stars on quasi-circular orbits. 

Fig. ^ compares the Plummer DF and the fitted DF 
and shows the estimated la deviation between both DFs. 


4 CONCLUSION 


As a conclusion we can state that the method can success¬ 
fully perform the iuversion of equation (^, given spectra of 
high enough quality (central signal to noise ratio equal to 
about 50 or higher) and the correct potential. The errorbars 
that can be derived from the Hessian matrix have a statisti¬ 
cal meaning and give a good idea of the uncertainty on the 
fitted model. 

This feature is probably the most important advantage 
of the present method. Up to now questions about the relia¬ 
bility of the inferred dynamics or the adopted potential have 
remained largely unanswered. With the present method, 
there is a fink with the S/N of the data, the choice of the 
template spectra and the continuum subtraction, which may 
all have critical repercussions on the inferred models. More¬ 
over, some simple form of population synthesis may be at¬ 
tempted for the best constrained models and geometries. 

In recent years, a lot of researchers have addressed the 
important subject of constraining the potential of invisible 
matter (i.e. dark halo matter, central black holes), using 
kinematical data (e.g. Merritt & Oh (tl996|) (who use ve¬ 
locity dispersion data), Rix et al. ( |l997D (who use and 
hi in addition to projected stream ing velocities and disper¬ 
sions) and Gerhard et al. (1997) (who use hi and veloc¬ 


ity dispersions)). It is not surprising that the inclusion of 
LOSVD shape parameters such as h^ and hi in a modeling 
procedure, in addition to projected streaming velocities and 
dispersions, gives a significant improvement in determining 
the potential. Since in the present method one uses in fact 
all the information on the shape of the LOSVDs that can 
be learned from the spectra, it should be a valuable tool to 
constrain the shape and extent of the potential. 


APPENDIX A: THE ANISOTROPIC 
VELOCITY MOMENTS 

The anisotropic velocity moments 
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Figure 8. The average LOSVDs for rp =0", 1", b" and 10^^ Experimental (black) and theoretical (grey) 1 cr errorbars are superposed. 
Where the errorbars coincide, the grey ones hide the black ones. Three positivity constraints are active. 


At2n,2m(j') = 27r / / F{E,L)Vr"v'^ dVrdVT, 


(Al) 


with L = rvT and E = tplr) — + vf ) where vt = 

are connected with the true velocity moments 

U2n,2i,2j{r) = / // F{E,L)vf^v^vl^ dvrdv4>dve (A2) 


(A3) 


through the relation 

U2n,2i,2j{r) = + ^)U2n,2(i+j)(r). 

It can be shown (Dejonghe 1986) that 

2"+"*r(n+i) 

fd'2n,2my) — -ttt-T ^ 


r(n + m) 


ri’i'r) 


{ip{r) — ^p{ip', r)) dip' (A4) 

do 

with p{ip, r) the augmented mass density, i.e. the density 
expressed as a function of r and ip, the potential. Inserting 


(P2|) for the mass density one hnds 


P>2n,2mlD — 


-ip(r 


DZ 




n+/3 


r(n + |)r(Q + i) 
r(a + n + m + 1) 


(1 + xY+'^ 


(AS) 


with X = (r jcPd'. Since 


DZ 


m+/3 


(1 + a:)'3+7 


= (-ir 


(1 +2;)'3+t' 


E 

3=0 


(—m — 




2/3 


2Fi(-m, 7 + /3; 1 + /3; 


r(m + /3+l) (r/c) 

r(/3 + l) (1 + (r/c)2)/3+7 

{rjcf , 


1 + {rjcf' 


(A6) 
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Figure 9. The same as in Fig. y, but now for a tangential q = —6 Plummer galaxy. The Fricke components with o = 7, 9, 10, 11, 12, 
13 and 14 are used in a fit to data with S/N ^ 80 for Vp = 0". 




Figure 10. The same as in Fig. but now for data from 
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tangential q = —6 Plummer galaxy with maximum S/N 


80. 
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Figure 11. The same as Fig. U but for a tangential q = —6 Plummer galaxy. 


the expression for the velocity moments can be rearranged 
as follows 

2 ^+^ r(n + i)r(a + l)r(/3 + m + 1 ) 

f^2n,2my) — - - - T7TT. 7TTTT- ^^^— X 




0? r(/3 + l)r(a + n + m + 1) 

iv 

\Q: + n+m \' / W ^ 


(1 + {r/cYY+~i 

{rtcf 


2 Fi(-m,7 + /3; 1 + /3; 


APPENDIX B: THE LOSVDS 


1 + (r/c)2 


(A7) 


To calculate the LOSVDs we make use of the fact that the 
Laplace transformation of a function j(x) can be written 
as a power series containing the moments of that function 

f*" = ir ■ 






(Bl) 


The aim is to calculate lp{r, rp, Vp), the probability of finding 
a star with projected velocity Vp at a distance r from the 
center of the cluster and on a line-of-sight with distance rp 
to the center of the galaxy, which has only non-zero even 
moments 


U 2 n{r,rp) = p{r){Vp{r,rp)} 


sr-^ ( 2n 


^ \ 2i 

i=0 ^ 
U 2 (n-i), 2 i,o(r} 


(cos (sin 77 )^’ X 


(B2) 


with rj(r, Vp) the angle between the spherical radial direction 
and the line-of-sight at a point of the line-of-sight. Using 
the above derived expressions for the velocity moments this 
becomes 


U2n{r,rp) = 


2/3 


2 ” r(n + |)r(Q-hi) (r/c) 

0r r(n-l-Q-l-l) {1 + {r / cY)!^+'^ 


] (-ir(sin77)^“ 

a—n / 

n 

E- 


z=0 

(-/3)i-i(/? + 7), 




j'Y-j)'- 


(B3) 


with ^ = (r/c)^/(l -I- (r/c)^). The Laplace transformation of 
lp(r, rp,Vp), for which we use the notation C, can be written 


r ^ r(n-f|) {-2'tpsY 

r(n-b a-b 1) n! \ i 

n=0 i=0 


“ rY-jy- 

00 j 00 . 

= Q^(/3 + 7)A(sinr7)^^X(-ir+^™ 

j^Q i=0 

-(-2V>s)*+^' X 


(sinr;)^" r(»+i + i) 


(i-bi)! r(i-bj-ba-bl) 

iFi{i-y j + i;i -bi -b a + 1; -2%j}s). 


(B4) 


with 


(B5) 


0F ^ (1 -b (r/c)2)/3+7 

Making use of the fact that 
Cyft{iFi{a-,b-, -y)} 

0<t<l (B6) 


r(a)r(&- a) 


= 0 


t > 1 


(B7) 


and 


C-Uy^Ct^y{f}} = + 
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Figure 12. The same as in Fig. g, but now for data from a radial q =1 Plummer galaxy with maximum S/N ^ 80. We used Fricke 
components with a = 4, 5, 6, 7 and 8. There are no active constraints. 





Figure 13. The same as in Fig. but now for a radial q = 1 Plummer galaxy. There are no active constraints. 
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Figure 14. The same as Fig. ^ but for a radial ^ = 1 Plummer galaxy. There are no active constraints. The data used in the fit have 
maximum S/N 80. 


k 

(B8) 

n=l 

where stands for the fc-th order derivative, one can in¬ 
vert the Laplace transformation. The function lp{r,rp,Vp) 
can be found by making the connection 

lp(r,rp,Vp) \vp\lp{r,rp,vl). (B9) 


The final result is 


lp{r,rp,Vp) = 


2/3 


1 r(Q-l-l) (r/c)_ 

y/^r{a + ^) (1-b (r/c)2)/3+7 




OO / 

T. I 


2 y (i + j)\ 

/ i+j '• 


(sin?7)^ 


„,2 \ -{i+j)+a-l/2 

1 _ 

2 tl>(r) 

1 „,2 


■2T’i(-(i-|-j),a; 


2 ’ 2ii{r 


(BIO) 


This formal result contains two infinite sums that diverge 
when implemented. Only when /3 is a positive integer and 
/3 -b 7 is a negative integer do both sums terminate after a 
finite number of terms. 


APPENDIX C: ERROR ANALYSIS ON THE 
COEFFICIENTS 


space, where the DF may reach its lower boundary b”, which 
is usually zero. Standard theory leads to the m so-called 
normal equations 

Dc -b AM = p (C2) 


that have to be solved, along with the r active constraints 
Ac = b, (C3) 

for the m-br unknowns c and d. One obtains the system 



This can be solved for the coefficients a 


a = M-/p^ + M-fb‘ 

= MDgignw" + Mpfb” 

= FPg^ + M-fh\ (C5) 

Here i and j range from 1 to m and s varies between 1 and 
r. The coefficients are written as a linear combination of the 
data points with the aid of the matrix PP = M~f gi^ (no 
summation over n), where n = 1,... ,N, which is indepen¬ 
dent of the data points gn. The covariance matrix of the 
coefficients, Sij(c), can be found using 

= {{Ci - Ci)iCj - Cj)} 

= FrFr{{9n-gn){gm-g^)) 

= ErE„^(3)F™. (C6) 


The minimisation problem is a least squares problem with 
constraints. The r active constraints can be dealt with by 
the introduction of r Lagrange multiplicators that form a 
vector d. The function that has to be minimised is 

= c*Dc — 2p*c -b e -b 2d*(Ac — b). (Cl) 

Here the element Ais of the matrix A is the value of the 
DF of the i-th component on the s-th grid-point in phase 


The covariance matrix of the data points is approximated 
by 'Sn-mig) = bnmI'uF. Inserting the explicit form of the 
elements of F one finds 

= (C7) 

In the absence of positivity constraints on the DF (A = 0) 
the covariance matrix of the coefficients is simply the inverse 
of the Hessian matrix. Taking the constraints into account 
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is expected to decrease the width of the errorbars since the 
positivity constraint rejects a lot of otherwise acceptable 
(but unphysical) fits to the data. 

Knowing the error on the coefficients it is also possible 
to obtain errorbars on every function that depends on these 
coefficients 

a(/(ci,...,c™))«ygE,,(c)^. (C8) 

For instance, the error on a LOSVD of the fitted model 
becomes 

G{(j){vp,Xp,yp)) = 

yjXpf yp'jJ^ij {'^pi Vp)- 

The elements of the Hessian matrix scale approximately pro¬ 
portional to the square of the signal to noise ratio. This 
means that the error on the coefficients, neglecting the in¬ 
fluence of the positivity constraints, scales inversely propor¬ 
tional to the signal to noise ratio. 

The Hessian matrix is also influenced by the fact 
whether or not the components are linearly independent. 
When two components are linearly dependent, their co¬ 
efficients can have an infinity of values without changing 
the outcome of the fit, balancing the one component with 
the other. The errorbars on the derived quantities, e.g. the 
LOSVDs, do not take into account this balancing and just 
measure the extreme values of these quantities, given the ex¬ 
treme values of the coefficients. Consequently, the errorbars 
will be infinitely large when linearly dependent components 
are used. This is reflected by the Hessian matrix being singu¬ 
lar, rendering its inversion impossible. When nearly linearly 
dependent components are used, the Hessian will be close 
to singular, resulting in large errorbars. 
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